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Abstract. An Ensemble Kalman Filter (EnKF, the predictor) is used 
make a large change in the state, followed by a Particle Filer (PF, the 
corrector), which assigns importance weights to describe a non-Gaussian 
distribution. The importance weights are obtained by nonparametric 
density estimation. It is demonstrated on several numerical examples 
that the new predictor-corrector filter combines the advantages of the 
EnKF and the PF and that it is suitable for high dimensional states 
which are discretizations of solutions of partial differential equations. 
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1 Introduction 

Dynamic Data Driven Application Systems (DDDAS) [1] aim to integrate data 
acquisition, modeling, and measurement steering into one dynamic system. Data 
assimilation is a statistical technique to modify model state in response to data 
and an important component of the DDDAS approach. Models arc generally 
discretizations of partial differential equations and they may have easily millions 
of degrees of freedom. The model equations themselves are posed in functional 
spaces, which arc infinitely dimensional. Because of nonlinearities, the probabil- 
ity distribution of the state is usually non-Gaussian. 

A number of methods for data assimilation exist [2]. Filters attempt to find 
the best estimate from the model state and the data up to the present. We 
present a combination of the Ensemble Kalman Filter (EnKF) [3] and the Se- 
quential Importal Sampling (SIS) particle filter (PF) [4]. The EnKF is a Monte- 
Carlo implementation of the Kalman Filter (KF). The KF is an exact method 
for Gaussian distributions. However, it needs to maintain the state covariance 
matrix, which is not possible for large state dimension. The EnKF and its vari- 
ants [6, 7] replace the covariance by the sample covariance computed from an 
ensemble of simulations. Each ensemble member is advanced in time by the 
model independently until analysis time, when the data is injected, resulting 



in changes in the states of the ensemble members. Particle filters also evolve a 
ensemble of simulations, but they assign to each ensemble member a weight and 
the analysis step updates the weights. 

The KF and the EnKF represent the probability distributions by the mean 
and the covariance, and so they assume that the distributions are Gaussian. 
This shows in the tendency of EnKF to smear distributions towards unimodal, 
as illustrated in Sec. 3 below. So, while the EnKF has the advantage that it 
can make large charges in the state and the ensemble can represent an arbitrary 
distribution, the EnKF is still essentialy limited to Gaussian distributions. On 
the other hand, the PF can represent non-Gaussian distributions faithfully, but it 
only updates the weights and cannot move ensemble members in the state space. 
Thus a method that combines the advantages of both without the disadvantages 
of either is of interest. The design of more efficient non-Gaussian filters for large- 
scale problems has been the subject of significant interest, often using Gaussian 
mixtures and related approaches [8]. 

The predictor-corrector filter presented here uses an EnKF as a predictor to 
move the state distribution towards the correct region and then a PF as corrector 
to adjust for a non-Gaussian character of the distribution. Nonparametric density 
estimation is used to compute the weights in the PF. The combined predictor- 
corrector method appears to work well on problems where either EnKF of PF 
fails, and it does not degrade the performace of the EnKF for Gaussian distribu- 
tions. Predictor-corrector filters were first formulated in [9, 10]. Related results 
and some probabilistic background can be found in [11]. 



2 Formulation of the Method 

A common procedure to construct an initial ensemble is as a sum with random 
coefficients [12], 

m 

u = ^ Xndnifn, dn N (0, 1) , {d„} independent, (1) 

n=l 

where {<p„} is an orthonormal basis in the space state V — M™ equipped with the 
Euclidean norm || • || . The elements of V are column vectors of values of functions 
on a mesh in the spatial domain. The basis functions (p„ are smooth for small 
n and more oscillatory for large n. If tlw^ tlw^ coefficients A„ sufficiently 
fast, the series (1) converges and u is a random smooth function in the limit as 
m — > oo. The sum (1) defines a Gaussian random variable with the eigenvalues 
of its covariance matrix equal to A^,. Possible choices of {^pk} include a Fourier 
basis, such as the sine or cosine functions, or bred vectors [2]. On the state space 
V, we define another norm by 

= U^^CniPn- (2) 

n=l n— 1 



Note that if = 1, \\-\\^ is just the original norm ||-|| on V. We generally use 
adapted to the smoothness of the functions in the initial ensemble, Xnji^n ~* 
as n — > oo. 

A weighted ensemble of N simulations {uk,Wk)^^i is initialized according 
to (1), with equal weights Wk = 1/N. The ensemble members are advanced by 
the model and at given points in time, new data is injected by an analysis step. 
The data consists of vector d of measurements, observation function h (u) — Hu, 
also called forward operator, here assumed to be linear, which links the model 
state space with the data space, and data error distribution, here assumed to be 
Gaussian with zero mean and known covariance R. The value of the observation 
function Hu is what the data vector would be in the absence of model and data 
errors. The value of the probability density of the data error distribution at 
the data vector d for a given value of the observation function Hu is called data 
likelihood and denoted by p{d\u). The probability distribution of the model state 
before the data is injected is called the prior or the forecast, and the distribution 
after the data is injected is called the posterior or the analysis. Assuming the 
forecast probability distribution has the density p^ , the density of the analysis 
is found from the Bayes theorem, 

p^iu) (xp{d\u)p^ (u), (3) 

where oc means proportional. 

Instead of working with densities, the probability distributions are approxi- 
mated by weighted ensembles. We will call the following analysis step algorithm 
EnKF-SIS. 

Predictor. Given a forecast ensemble 

JV ^ 



k=l 

the members of the analysis ensemble are found from the EnKF, 

ul = ui + K{dk-Hui), dk~N{d,R), K = QH^ {HQH^ + 

where dk are randomly sampled from the data distribution, and Q is the forecast 
ensemble covariance, 

N N 

Q = '^Wk {uk -u-l) [uk -u^) , ^^^wlu{. (4) 
fc=i fe=i 

This is the EnKF from [3], extended to weighted ensembles by the use of the 
weighted sample covariance (4). 

Corrector. The analysis members are thought of as a sample from some 
proposal distribution, with density p^ . Ideally, the analysis weights should be 
computed from the SIS update as [4] 



ENKF 




Fig. 1. Data assimilation witli bimodal prior. EnKF fails to capture the non-Gaussian 
features of the posterior, but both SIS and EnKF-SIS represent the nature of the 
posterior reasonably well. 



However, the ratio of the densities is not known, so it is replaced by a nopara- 
metric estimate inspired by [13], giving 

Ell f 11 1i/ N 

ui-u? <hk k -f-^ 

2^e:\\u1-ul\\^<hk N fc=i 

The bandwidth hk is the distance from to the [A^^/^J-th nearest member u", 
measured in the norm. 

3 Numerical Results 

Fig. 1 demonstrates a failure of EnKF for non-Gaussian distributions, while SIS 
and EnKF-SIS do fine. We construct a bimodal prior in ID by first sampling 
from N (0, 5) and assigning the weights by 




Fig. 2. Ensemble filters mean and optimal filter mean for stochastic ODE (5). EnKF- 
SIS was able to approximate the optimal solution better than either SIS or EnKF 
alone. 



The data likelihood is Gaussian. The ensemble size was A'' = 100. 

The next ID problem demonstrates that EnKF-SIS is doing better that either 
EnKF or SIS alone in filtering for the stochastic differential equation [14] 

— = 4u - 4u'^ + KT], (5) 

where r]{t) is white noise. The parameter k controls the magnitude of the stochas- 
tic term. 

The deterministic part of this differential equation is of the form 

du , 

where the potential /(w) = —2u^ + u^. The equilibria are given by /' (w) = 0; 
there arc two stable equilibria at ?i = ±1 and an unstable equilibrium at u = 0. 
The stochastic term of the differential equation makes it possible for the state to 
flip from one stable equilibrium point to another; however, a sufficiently small k 
insures that such an event is rare. 

A suitable test for an ensemble filter will be to determine if it can properly 
track the model as it transitions from one stable fixed point to the other. From 
Fig. 1, it is clear that EnKF will not be capable of describing the bimodal nature 
of the state distribution so it will not perform well when tracking the transition. 
Also, when the ensemble is centered around one stable point, it is unlikely that 
some ensemble members would be close to the other stable point. It is known 
that SIS can be very slow in tracking the transition and EnKF can do better 
[14]. Fig. 2 demonstrates that EnKF can outperform both. 



The solution u of (5) is a random variable dependent on time, with density 
p{t, u). The cvohition of the density in time is given by the Fokker-Planck equa- 
tion, which was solved numerically on a uniform mesh from u = — 3 to u = 3 
with the step Au = 0.01. At the analysis time, the optimal posterior density 
was computed by multiplying the probability density of u by the data likelihood 
following (3) and then scaling so that again / pdu = 1, using numerical quadra- 
ture by the trapezoidal rule. The data points were taken from one solution of 
this model, called a reference solution, wliic;li exhibits a switch at time t « 1.3. 
The data error distribution was normal with the variance taken to be 0.1 at 
each point. To advance the ensemble members and the reference solution, we 
have solved (5) by the explicit Euler method with a random perturbation from 

1 /2 

N{0,{At) ' ) added to the right hand side in every step [16]. The simulation 
was run for each method with ensemble size 100, and assimilation performed for 
each data point. 

Finally, typical results for filtering in the space of functions on [0, tt] of the 
form 

d 

u = c„ sin {nx) (6) 

arc in Figs. 3 and 4. The ensemble size was TV = 100 and the dimension of 
the state space was d = 500. The Fourier coefficients were chosen A„ = to 
generate the initial ensemble from (1), and «;„ = for the norm in the density 
estimation. 

Fig. 3 again shows that EnKF cannot handle bimodal distribution. The prior 
was constructed by assimilating the data likelihood 

C 1/2 if u{-k/A) and 
p{d\u)=l M(37r/4)e (-2,-1) U (1,2) 
I otherwise 

into a large initial ensemble (size 50000) and resampling to the obtain the forecast 
ensemble size N = 100 with a non-Gaussian density. Then the data likelihood 
u {tt/2) — 0.1 ~ N{0, 1) was assimilated to obtain the analysis ensemble. 

Fig. 4 shows a failure of SIS. The prior ensemble sampled from Gaussian 
distribution with coefficients A„ = using (1) and (6), and the data likelihood 
was u(7r/2) -7~ Ar(0,l). 

4 Conclusion 

We have demonstrated the potential of a predictor-corrector filter to perform a 
successful Bayesian update in the presence of non-Gaussian distributions, large 
number of degrees of freedom, and large change of the state distribution. Open 
questions include convergence of the filter in high dimension when applied to 
multiple updates over time, mathematical convergence proofs for the density 
estimation and for the Bayesian update, and performance of the filters when 
applied to systems with a large number of different physical variables and modes, 
as is common in atmospheric models. 



Prior Ensemble 



Posterior: EnKF 




Fig. 3. EnKF smears non-Gaussian distribution. The horizontal axis is the spatial co- 
ordinate X £ [0, 7r]. The vertical axis is the value of function u. The level of shading on 
each vertical line is the marginal density of -u at a fixed x, computed from a histogram 
with 50 bins. While EnKF completely ignores the non-Gaussian character of the pos- 
terior and centers the distribution around u = 0, EnKF-SIS shows darker bands at the 
edges. 
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